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1.  Introduction 


In  this  report  we  consider  the  problem  of  finding  all  the  zeros  of 
the  polynomial 

f(z)  =  a^  +  a^z  +...+a^z"  (1.1) 

and  estimating  error  bounds  for  them.  To  find  the  set  of  zeros  we  use 

the  algorithm  of  Madsen  (1973),  which  we  have  found  to  compare  favourab 

/ 

with  other  algorithms.  For  error  estimation  we  apply  Rouche's  theorem 
as  recommended  by  Peters  and  Wilkinson  (1971),  but  with  some  difference 
in  detail.  Both  of  these  algorithms  are  a  little  simpler  in  the  case 
where  the  polynomial  has  complex  coefficients,  so  we  describe  the 
algorithms  for  this  case  in  sections  2  and  3  and  the  modifications  for  i 
real  case  in  section  4. 

We  believe  that  our  code  is  in  accord  with  the  ANSI  standard,  havir 
checked  it  with  Bell  Telephone  Laboratories' Fortran  verifier  (Ryder, 
1973)  and  run  our  test  programs  with  array  subscript  checking.  We 
describe  this  code  in  section  5  and  in  an  appendix  give  specification 
sheets,  and  listings  (produced  by  the  Bell  Laboratories' verifier,  so  th= 
include  cross  references).  For  the  Harwell  subroutine  library  we  have 
made  a  small  number  of  changes  in  order  to  shorten  the  argument  lists, 

t 

at  the  expense  of  a  departure  from  ANSI  standard.  We  have  decided  aga: 
including  a  single-length  version  in  the  library  because  the  IBM  370/16 
has  a  very  short  single-length  word  (6  hexadecimal  digits,  so  that 
numbers  just  greater  than  unity  are  held  to  about  1  part  in  10^)  but  he 
double-length  hardware  which  executes  very  little  slower  than  the  singl 
length  hardware.  Some  results  obtained  with  the  library  subroutines, 
together  with  comparisons  with  other  algorithms  are  given  in  section  6. 
Our  code  (listed  in  section  7)  contains  comments  suitable  for  machine 
processing  which  detail  the  changes  needed  for  the  single-length  and 


Harwell  li.'rary  versions. 


We  would  like  to  acknowledge  the  help  of  M.J.D.  Powell  in  carefully 
checking  a  draft  of  ;;his  report  and  making  several  valuable  suggestions. 

2.  Finding  the  routs  iti  the  complex  case 

We  use  the  algoncnm  of  Madsen  (1973)  to  find  the  .'oot  of  minimal 
or  near-minimal  modulus  and  then  use  forward  deflation  to  construct  a 
polynomial  of  degree  n-1  whose  roots  are  the  remaining  roots  of  the 
original  polynomial.  The  process  is  repeated  until  approximations  to 
all  the  roots  have  been  found.  Wilkinson  (1963)  has  shown  that  forward 
deflation  is  stable  provided  a  large  root  is  not  accepted  before  a  much 
smaller  one.  Our  algorithm  does  not  guarantee  that  the  moduli  of  the 
roots  are  strictly  increasing  but  our  experience  has  always  been  that  they 
are  found  in  roughly  increasing  order.  A  version  with  the  composite 
deflation  of  Peters  and  Wilkinson  (1971),  which  should  be  stable  no  matter 
in  which  order  the  zeros  are  found,  was  tried  but  it  did  not  give  more 
accurate  results.  In  any  case  it  is  not  clear  how  to  apply  the 
composite  deflation  when  two  complex  conjugate  roots  of  a  real  polynomial 
have  been  found  so  that  deflation  by  a  real  quadratic  factor  is  wanted. 

It  remains  necessary  to  describe  Madsen's  (1973)  algorithm  in  detail 
and  we  will  consider  its  application  to  the  original  polynomial .  The 


general  strategy  of  the  algorithm  is  that,  given  an  iterate  Z|^,  a 
tentative  step  dZj^  is  found  and  the  next  iterate  Z|^^.|  is  taken  at  the  best 
point  (in  the  sense  of  |f(z)|)  encountered  in  a  short  search  of  values 


on  the  line  through  Zj^  and  Zj^+dZj^.  Because  the  search  may  sometimes 
yield  no  better  value  than  that  at  Z|^  we  may  sometimes  have  Z|^^.|=Z|^  and 
in  this  case  ensure  that  the  next  tentative  step  is  shorter  and  in  a 


different  direction.  The  inclusion  of  searches  ensures  rapid  convergence 


to  multiple  roots  and  reliable  convergence  when  difficulties  are 


ATTRIBUTES:  Col  1 :  C  if  in  COMMON 


encountered.  Such  a  search  is,  however,  wasteful  if  we  are  so  near  a  simple 
root  that  Newton's  iteration  is  reliable  and  fast.  We  have  therefore 

devised  a  test  (given  by  inequality  (2.7))  which  normally  ensures  this. 

While  the  algorithm  performs  searches  we  say  it  is  in  stage  1;  otherwise 
it  is  in  stage  2,  performing  straightforward  Newton  iteration.  It  begins  in 
stage  1,  which  we  now  describe. 

The  tentative  step  dZj^  is  found  with  the  help  of  stored  values  of 
^k’  ^'^^k^’^k-1  previous  tentative  step  d2|^_i-  If  the 

last  iteration  was  successful  (Z|^  ^  ^k-1  ^  Newton  correction 

=  -f(z^)/f'{z^)  (2.2) 

is  calculated  and  the  next  tentative  step  is  taken  as 

l"kl  "  3|2,^-z^.^|  (2.3a) 

[_3|Z|^-Z|^_^  |e’®  n^/lnj^l  otherwise  (2.3b) 

where  6  is  chosen  (rather  arbitrarily)  as  arctan(3/4).  If  the  last 
step  was  unsuccessful  (Zj^*Z|^_l)  then  we  take  the  tentative  step  to  be 

dz^  =  -J  e^®  dz,^.^.  (2.3c) 

After  a  successful  iteration  we  normally  expect  to  want  to  take  a  Newton 
step  (2.3a),  but  we  include  the  alternative  (2.3b)  because  accidently 
coming  near  to  a  stationary  point  of  f(z)  is  likely  to  make  the  Newton  step 
ridiculously  large.  We  include  a  change  of  direction  in  (2.3b)  because 
if  a  saddle  point  is  being  approached  the  direction  nj^  may  be  a  worse 
search  direction  than  almost  any  other.  After  an  unsuccessful  iteration 
we  want  to  change  the  search  direction  to  one  likely  to  be  successful 
and  reduce  the  step  size;  this  leads  to  formula  (2.3c).  Its  repeated 
use  is  sure  to  yield  a  descent  direction. 
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The  algorithm  used  by  PA06BD  to  find  the  roots  has  already  been 


Once  the  tentative  step  has  been  found  we  test  the  inequality 


|f(z,^+dz,^)l  <  |f(Z|^)l  .  (2.4) 

If  this  is  satisfied  then  we  calculate  the  numbers 

lf(Z|^+p  dZ|^)(,  p=l,2,...,n  (2.5) 

continuing  for  as  long  as  these  are  strictly  decreasing.  If  inequality 
(2.4)  does  not  hold  th  .n  we  calculate  the  numbers 

|f(z,^  +  dz^/2P)l.  p=0.1,2.  lf(z,^  +  ie^®  dz^)l  (2.6) 

again  continuing  until  the  sequence  ceases  to  decrease.  In  all  cases  we 
take  Z|^^i  to  be  the  best  point  found.  Note  that  if  there  is  a  true  multiple 
zero  of  multiplicity  m  or  if  we  are  a  fair  distance  from  a  cluster  of  m 
zeros  then  Z|^+mn|^  will  be  a  very  good  estimate  of  the  solution  and  will  be 
found  by  our  search.  In  fact  we  get  quadratic  convergence  to  a  multiple 
zero.  Note  also  that  when  an  iteration  fails  following  a  search  to  the 
end  of  sequence  (2.6),  the  choice  (2.3c)  leads  to  a  step  likely  to  be  in 
a  direction  of  decreasing  )f(z)l. 

To  complete  our  description  of  stage  1  we  need  to  specify  starting 
iterates.  We  take  these  to  be 


z 


0 


=  0 


^1 


_f-f(0)/f’(0)  if  f’(0)=0 
~l  1  otherwise 

dz. 

J  mir 


IdZol 


(2.7) 


The  iteration  really  starts  from  z^  but  we  have  to  include  z^  and  dz^ 
because  they  are  needed  for  choosing  the  tentative  step  dz-j .  This  choice 
of  z^  is  used  because  its  modulus  is  certainly  less  than  that  of  any  root 
of  f(z)  and  it  is  in  the  direction  of  steepest  descent  of  |f(z)|  from  the 
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origin.  It  is  therefore  likely  that  we  will  converge  to  a  root  of  near- 
minimal  modulus. 


We  do  not  make  our  main  test  for  switching  to  stage  2  (straightforward 
Newton  iteration)  until  a  stage  1  search  has  led  to  the  choice 
Zk^l=Zj^+d2k.  This  is  obviously  sensible  and  has  the  added  virtue  that 
we  should  (correctly)  avoid  the  switch  when  converging  to  a  multiple  root. 
The  test  itself  is  based  on  the  Kantorowitz  theorem  (see,  for  example, 
Ostrowski  (1966))  which  states  that  if  is  the  circle  with  centre 
Z|^+n|^  and  radius  Inj^j  where  nj^  is  the  Newton  step  (2.2)  then  the 
conditions 


f(z^)  f(z^)  0 

2lf(z^)|  max  |f"(z)|  <  \r 


ensure  the  convergence  of  Newton's  iteration  starting  from  Zj^.  This 
leads  us  to  test  the  inequality 


2|f(z^)|lf'(z^.^)  -  f'(z^)|  s  |f'(z,^)r  Uk-r^kl 


Of  course  this  is  not  equivalent  to  test  (2.8)  because  we  have  replaced 

max  |f"(z)|  by  a  rather  crude  difference  approximation  but  we 
zeK 

0 

nevertheless  expect  it  usually  to  predict  correctly  that  straightforward 
Newton  iteration  will  be  satisfactory.  We  check  inequality  (2.9)  at  every 
step  in  stage  2  and  switch  back  to  stage  1  if  it  is  violated.  We  also 
check  the  inequality  (2.4)  and  if  this  is  violated  return  to  stage  1 
beginning  by  modifying  the  tentative  step  with  formula  (2.3c)  as  in  stage  1. 

We  complete  this  section  by  describing  our  convergence  criterion. 

We  terminate  if  a  stage  1  search  or  a  stage  2  iteration  leads  to  a  new 
iterate  Zj^^i  different  from  Zj^  and  yet  such  that  the  inequality 


i^+r"kl  " 

holds  where  e  is  the  largest  number  such  that  to  machine  accuracy  1+c  =  1. 

We  also  terminate  if  the  condition 

l^(Zk+l)l  =  -  ISnlaJe  (2.11) 

holds.  The  expression  Ibnla^le  is  a  generous  overestimate  of  the  final 
roundoff  made  in  calculating  f(z)  at  the  root  of  smallest  modulus  and  we 
expect  that  such  accuracy  will  be  attainable.  The  normal  convergence 
pattern  is  that  |f(Z|^)|  decreases  until  well  below  lenja^le  and  then 
roundoff  errors  cause  a  new  iterate  Z|^^^=Z|,  to  be  taken  so  that  (2.11)  is 
satisfied.  If  such  accuracy  is  unattainable  then  the  step  will  decrease 
steadily  because  of  the  application  of  (2.3c)  until  (2.10)  is  satisfied. 

This  combination  of  convergence  criteria  means  that  we  are  certain  to 
obtain  a  good  solution  and  almost  certain  to  obtain  the  best  possible. 
Furthermore  this  result  is  usually  obtained  with  only  one  more  iteration  than 
is  necessary  to  get  this  good  accuracy. 

3.  Error  estimation 

We  seek  to  estimate  error  bounds  for  all  the  roots  produced  by  the 
algorithm  of  the  previous  section.  We  base  our  algorithm  on  'tiir'.t  cf 
Peters  and  Wilkinson  (1971)  and  will  follow  their  notation.  The  most 
significant  difference  is  that  they  look  for  non-overlapping  discs  each  of 
which  contains  precisely  the  same  number  of  exact  and  approximate  roots, 
whereas  we  allow  overlapping  discs.  This  's  because  it  can  happen  that 
one  root  is  well  determined  although  it  is  inside  the  best  disc 
obtainable  for  another  (ill -determined)  root.  We  therefore  look  for  a 
separate  disc  for  each  root  and  take  its  centre  to  be  at  the  calculated  root 
itself.  This  also  allows  a  very  simple  form  of  output  to  the  user  since 
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all  that  is  required  is  a  radius  for  each  root. 

We  suppose  that  we  have  approximate  roots  a.,  i=l,2,...,n.  Then 
the  polynomial 

P(z)  =  a  n  (z-a.)  (3.1) 

"  i=l  ^ 

should  agree  with  the  original  polynomial  (1.1),  but  because  the  roots 
are  not  exact  there  will  be  an  error 

Q(z)  =  P(z)  -  f(z).  (3.2) 

Now  Rouche's  theorem  states  that  if  P(z)  and  Q(z)  are  analytic 
functions  in  and  on  the  closed  curve  C  and  the  inequality 

|Q(z)I  <  tP(z)l  (3.3) 

holds  on  C  then  P(z)  and  P(z)-Q(z)  have  the  same  number  of  zeros  inside  C. 
We  apply  this  here  by  looking  ^or  circles  with  centres  ,  i=l,2,...,n  on 
which  condition  (3.3)  holds.  The  following  theorem  shows  that  there  is 
no  need  to  worry  about  the  overlapping  of  some  of  these  discs. 

Theorem  1  If  condition  (3.3)  holds  on  each  of  the  circles  centre  , 
radius  r^.,  i=l,2,...,n,  then  the  roots  (p.  of  f(z)  may  be  ordered  so 
that 

5  r^.,  i=l,2,...,n  .  (3.4) 

Proof  Regard  the  perimeters  of  the  circles  as  dividing  the  plane  into  a 
set  of  non-overlapping  regions  R^,  each  of  which  is  the  intersection  of  a 
subset  of  the  set  of  discs  that  the  circles  enclose  and  their  complements. 
A  simple  case  is  illustrated  in  Figure  1.  Let  R.  be  the  region 
containing  a^,  i=l,2,...,n.  Note  that  a  region  may  contain  more  than  one 
a^.  so  that  there  may  be  coincidences  among  the  (e.g.  ki=k2=3  in 
Figure  1).  The  set  of  regions  R.  ,  i=l,2,...,n  together  contain  all  the 


'1  R 


Figure  1  Four  circles  centre  a^.  dividing  plane  into  seven  regions  R. 


roo^s  a.  of  P(z)  and  each  R  has  a  boundary  on  which  condition  (3.3) 
holds  so  contains  exactly  the  same  number  of  roots  of  f(z)  as  of  P(z). 
Therefore  the  regions  R  (i=l contain  all  the  roots  of  f(z)  and 
these  may  be  ordered  so  that  R.  contains  ({).  (i=l  ,2,...,n).  The  result 

Ki  1 

now  follows  since  each  R,  is  contained  in  the  disc  centre  a.  radius  r-. 

Being  able  to  use  overlapping  discs  leads  to  simplifications  in 
coding  and  sometimes  to  much  better  error  bounds.  This  improvement  is 
illustrated  by  the  example  shown  in  Figure  1,  where  and  are  quite 
ill-conditioned  but  O4  ^re  not.  The  procedure  of  Peters  and 

Wilkinson  would  have  forced  us  to  regard  all  four  a.  as  a  cluster  and  to 
enclose  them  in  a  single  disc.  This  would  have  made  little  difference  to 
the  bounds  for  and  02  but  would  have  significantly  worsened  the  bounds 
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for  and  a^.  Another  illustration  of  this  point  is  given  in  §6  with 
an  example  using  our  code. 

We  now  suppose  that  the  error  polynomial  is 


Q(z)  =  I 

1=0 


(3.5) 


and  describe  how  to  find  a  suitable  circle  for  a  typical  root.  For 
simplicity  of  notation  let  us  imagine  that  the  roots  are  reordered  so 
that  the  one  for  which  we  are  seeking  a  circle  is  and  |a^.-a^[  increases 

monotonically  with  i.  Rouche's  condition  (3.3)  is  satisfied  on  a  circle 
centre  a-j  and  radius  r  if  the  inequality 

n  .  n 


I  leJ  (r+laJ)^  <  |aj,  n  (  |  -  r)|  (3.6) 

i=o  i=l 


holds.  If  m  is  such  that 


iv“ii  ^  ^  ivr“i 


(3.7) 


(we  ignore  the  case  r  =  |ot.-a^|  for  some  i  because  clearly  inequality  (3.6) 
cannot  be  satisfied  in  this  case)  then  we  may  rewrite  inequality  (3.6)  in 
the  form 


k(r)  <  (r  -  |) 


(3.8) 


where  k(r)  is  given  by  the  equation 

j  (r+lct^l)^’ 


k(r)  = 


l^nl  JT  (l“i-“il-0 

i=m+l 


(3.9) 


Peters  and  Wilkinson  (1971)  solve  the  equation 


(1.1)  k(0)  =  (r,  -  max  1)" 

l<i<m 


(3.10) 
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and  then  check  whether  r-j  satisfies  inequalities  (3.7)  and  (3.8).  Almost 
always  inequality  (3.8)  holds  because 

(i)  the  right  hand  side  of  (3.10)  underestimates  the  corresponding 
expression  in  (3.8), 

(ii)  the  safety  factor  (1.1)  has  been  introduced 

and  (iii)  k(r)  is  nearly  constant  in  the  usual  case  where  the  roots 

i=l,2,...,m  are  well  separated  from  the  rest. 

They  give  no  recommendation  for  dealing  with  the  case  where  r-j  does  not 

satisfy  (3.8),  but  presumably  intend  that  an  iteration  should  be  set  up 

with  0  and  r^  in  (3.10)  being  replaced  by  r^  and  .  This  yields  a 

montonically  increasing  sequence,  so  that  we  must  eventually  either 

satisfy  inequality  (3.8)  or  break  the  right-hand  inequality  (3.7) 

necessitating  a  new  start  with  a  greater  value  for  m.  Actually  Peters  and 

m 

Wilkinson  take  the  centre  of  their  circle  to  be  a  =  Z  a-/m  rather  than 

i=l  ^ 

a.|  and  this  means  that  the  right-hand  side  of  (3.10)  is  quite  a  good 
estimate  of  the  right-hand  side  of  (3.8)  so  that  the  procedure  gives  a 
realistic  radius  r  without  very  much  computation.  Unfortunately  this  is 
not  the  case  with  a-j  in  use,  except  when  m=l  and  we  have  sometimes 
obtained  solutions  such  that  the  right-hand  side  of  (3.8)  exceeds  the  left 
by  factors  as  big  as  1,000.  We  have  therefore  decided  to  use  an  iteration 
based  directly  on  (3.8).  Given  an  iterate  r.  we  seek  a  new  iterate  such 

vl 

that 

m 

(1.05)  k(rj)  <  (r^^^  -  {a^.-a^D  <  (1.1)  k(rj).  (3.11) 

In  view  of  inequality  (3.7)  it  seems  sensible  to  start  with  r^  =  l«^“ail 
rather  than  Wilkinson  and  Peters'  r„=0.  To  solve  (3.11)  we  use  the  method 
of  bisection  because  m  is  usually  small  so  that  function  evaluations  are 
cheap,  no  great  accuracy  is  required  in  view  of  the  slack  in  (3.11), 


# 


and  suitable  initial  upper  and  lower  bounds  are  available  in  r^  and  the 
Peters  and  Wilkinson  overestimate  obtained  from  the  analogue  of  equation 
(3.10).  Typically  between  four  and  seven  iterates  are  required. 

Peters  and  Wilkinson  suggest  that  r^^^  should  be  checked  directly 
in  (3.8)  but  a  simple  test  often  avoids  any  need  for  this.  If  the 


inequality 


^j^l  "  1^1 
Tj  +  la,! 


'“m+1  '  *^1 1  ~ 

l“m+l  '  “it  - 


<  1.05 


(3.12) 


holds  then  it  follows  from  the  definition  (3.9)  that  the  inequality 


k(rj^^)  <  1.05  k(rj.)  (3.13) 

also  holds.  From  this  last  inequality  and  the  first  inequality  (3.11) 

we  deduce  that  rj^^  satisfies  inequality  (3.8). 

Rounding  errors  occur  when  the  coefficients  of  P(z)  are  calculated, 

but  can  be  reduced  by  multiplying  out  the  factors  (z-ct^. )  in  order  of 

increasing  |a^l.  We  therefore  use  the  same  order  as  that  in  which  they 

were  calculated.  If  bounds  e^  on  the  errors  were  available  we  could 

n  i 

work  with  the  error  polynomial  E  (leJ+e.)z  in  place  of  the  polynomial 

i=o  '  ^ 

(3.5)  and  obtain  strict  bounds  on  the  errors  in  the  calculated  roots. 

We  hoped  to  use  the  running  error  analysis  of  Peters  and  Wilkinson  for 
this  purpose,  but  unfortunately  found  that  it  sometimes  gave  gross  over¬ 
estimates,  as  explained  in  the  last  paragraph  of  this  section.  We 
therefore  have  ignored  this  source  of  error,  relying  on  such  errors  being 

reflected  in  larger  coefficients  e^.  in  the  polynomial  Q(z).  It  should  be 

r 

noted  that  the  sequence  of  constructed  polynomials  IT  (z-a-)»  r=l,2,...,n 

i=l  ^ 

is  not  identical  with  the  sequence  of  deflated  polynomials  used  when 
finding  the  roots  because  we  perform  the  multiplication  in  the  same  order 
as  that  in  which  the  roots  were  found.  Therefore  errors  produced  by  the 


multiplication  are  likely  to  show  in  enlarged  coefficients  e^.  Not 
bounding  these  errors  means  that  our  bounds  are  not  strict  bounds,  but 
they  are  more  realistic.  In  none  of  our  tests  did  the  actual  errors  exceed 
the  estimated  bounds. 

It  is  straightforward  to  allow  for  the  effects  of  uncertainties  in 
the  original  coefficients  by  working  with  the  error  polynomial 

n 

iio  b^,  i=0,l,...,n,  are  bounds  on  these  errors.  This 

we  do  in  our  subroutine. 

We  complete  this  section  by  explaining  why  the  running  error 

n 

analysis  of  Peters  and  Wilkinson  for  the  coefficients  of  n  (z-a.) 

i=l  ^ 

sometimes  gives  very  pessimistic  results.  If  the  computed  product 

n  (z-oj)  is  1  b('"^z’  then  they  generate  coefficients  fl*"^  from  the 
i=l  i=o  '  ' 

recurrences. 

fi""  =  *  iviifi'"’  *  iviif’S’'’!  '■■’•2 . <'•”> 


i=l  ,2,...,r 


(r+1)  .  (r+1)  _ 

-1  ■  V+1  "  ^ 


to  yield  bounds  2  for  the  coefficients  of  P(z)  if  floating-point 

computation  with  t  binary  places  is  used.  This  can  sometimes  lead  to 
very  pessimistic  bounds  since  the  recurrence  (3.14)  fails  to  distinguish 
properly  between  n(z-a-)  and  n(z+|a^.l).  In  fact  the  recurrence  (3.14) 
is  obtained  by  bounding  the  corresponding  recurrence 

^i  ®i-l  “r+1  ®i  ^l“r+ri  ^2  i  •  =  ; 


NPl  (standard  versions  only)  is  an  INTEGER  which  must  be  set  to  n+1 . 


Such  bounds  are  particularly  unrealistic  in  such  a  case  as 
f{z)  =  z^-1  since  the  numbers  grow  with  r  in  much  the  same  way  as  the 
recurrence  for  generating  the  coefficients  of  (z+l)"  and  these  coefficients 
are  '^Cr,r=0,l . ,n  .  In  fact  the  recurrence  (3.14)  has  some  added 
terms  when  compared  with  the  recurrence  for  the  coefficients  of  (z+1)”. 

4.  Changes  to  the  algorithms  in  the  real  case 

The  only  change  made  to  our  root-finding  algorithm  in  the  case  where 
f{z)  has  real  coefficients  is  that  once  a  complex  root  has  been  found  it  is 
either  perturbed  into  a  real  root  or  its  conjugate  is  taken  as  another  " 
This  ensures  that  each  deflated  polynomial  is  also  real  and  that  work  i 
saved  for  genuinely  complex  roots.  Once  a  complex  root  has  been 

found  we  evaluate  f(X|^)  and  use  Peters  and  Wilkinson's  (1971)  running 
error  analysis  to  bound  the  roundoff  error  made  in  this  evaluation.  The 

recurrence  used  for  evaluating  ^(Xj^)  is  given  by  the  equations 


*  a„ 
n  n 


^i'Vi+i-^^i’  . 0 

s. 


(4.1) 


-0 

and  the  corresponding  running  error  analysis  is  given  by  the  equations 


9n  =  0  1 

}  (4.2) 

9i  "  I’^kl^^i+i  +  ls^+ll)  +  Is-l,  i=n-l,...,0  J 
to  yield  the  bound  eg^  for  the  error  in  f(X|^)  where  e  is  the  relative 
floating-point  accuracy.  If  the  inequality 

hoi  <  2€g^  +  (4.3) 


holds  (where  f^(a(^)  is  the  computed  value  at  aj^),  we  take  v,  to  be  a  real 
root  and  deflate  with  it.  Otherwise  we  deflate  with  the  complex  conjugate 

pair  X|^  ±  iyj^.  It  is  very  important  that  the  case  where  we  have  a  simple 


real  root  well  separated  from  the  other  roots  should  not  be  taken  as  a 
complex  pair  because  a  double  deflation  in  such  a  case  would  be  disastrous. 

We  therefore  consider  this  case  and  suppose  that  ({)|^  is  such  a  real  root  and 
is  the  (complex)  approximation  to  it  found  by  the  algorithm.  Xj^  (real 
part  of  is  a  better  approximation  to  (t>|^  in  the  sense  that  the 
inequality 

holds  since  <})|^  is  real.  Further  the  exact  value  f(z)  behaves  like  a 
constant  times  (z-<{i|^)  if  z  is  near  It  is  therefore  reasonable  to  expect 

that  the  inequality 

lf(Xk)l  ^  (4.5) 

holds  between  the  corresponding  exact  values.  We  therefore  expect  Xj^  to  be 
at  least  as  good  a  root  as  so  that  it  is  virtually  certain  that  in¬ 
equality  (4.3)  will  hold,  since  eg^  is  a  rigorous  upper  bound  on  the  round¬ 
off  error  in  computing  f(X|^)  and  may  also  be  taken  (slightly  incorrectly) 
as  a  bound  on  the  error  in  computing  |f{ot(^)|. 

A  disadvantage  of  the  test  (4.3)  is  that  we  may  decide  we  have  a 
real  root  when  a  complex  conjugate  pair  is  the  true  position.  We  feel 

that  this  disadvantage  is  quite  slight  because  it  is  reasonable  to  regard 
any  point  for  which  inequality  (4.3)  holds  as  a  good  approximate  root. 

The  only  change  to  the  error  bounding  algorithm  is  that  once  a 
bound  for  a  complex  root  has  been  found,  this  bound  is  also  used  for  its 
complex  conjugate. 
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5 .  Descripti on  of  Fortran  subroutines 


0 

In  this  section  we  describe  the  Fortran  subroutines  themselves.  It  is 
a  short  section  because  we  have  included  rather  full  comments  in  the  code 
itself,  since  we  believe  that  this  provides  the  most  convenient  form  of 
documentation.  Specification  sheets  and  listings  of  the  code  are  given 
in  section  7.  The  specification  sheets  document  all  three  versions  (double- 
and  singe-length  standard  Fortran  and  double-length  IBM  Fortran).  The 
listings  are  of  the  standard  Fortran  double-length  versions  and  contain 
specially  coded  comments  detailing  changes  needed  for  the  other  versions. 

We  have  been  using  a  simple  preprocessor  (written  in  standard  Fortran)  to 
convert  from  one  version  to  another.  Where  a  statement  differs  between 
versions  we  preceded  the  genuine  (double-length  standard)  statement  by  the 
alternative  version  or  versions  modified  by  the  insertion  of  C  in  column  1, 

/  in  column  72,  and  the  letter  I  or  S  (for  IBM  version  or  single-length 
version)  in  column  71. 

The  listings  used  are  those  produced  by  the  Bell  Telephone 
Laboratories’ Fortran  verifier  because  they  include  cross-references  which  are 
very  helpful  when  reading  the  code.  They  consist  of  lists  of  all  the 
identifiers  and  labels  in  lexographical  order  together  with  the  statement 
numbers  of  all  references  to  them  and  the  following  coded  information 
for  each  identifier : 

TYPE:  Col  1:  E  if  explicitly  typed 

Col  2:  I  for  integer 
P>  for  real 

D  for  double  precision 
C  for  complex 
L  for  logical 
H  for  Hollerith 

USE:  FA  for  arithmetic  statement  function  argument 

FN  for  function 

E  for  external  subroutine  or  function 
GT  for  assigned  "go  to"  variable 
IF  for  intrinsic  function 
SN  for  subroutine 
V  for  variable 

-15- 


1 


ATTRIBUTES:  Col  1 :  C  if  in  COMMON 

Col  2:  E  if  an  EQUIVALENCE 
Col  3:  A  if  a  dummy  argument 
Col  4:  S  if  value  set  by  program  unit 
Col  5:  S  if  scalar 
A  if  array 

Col  6:  if  array  then  number  of  dimensions 

Because  double-length  complex  facilities  are  not  available  in 
standard  Fortran  we  do  not  use  any  complex  variables,  although  the  arguments 
A,R  and  E  may  be  regarded  as  C0MPLEX*16  (or  COMPLEX  for  the  single-length 
version)  on  the  IBM  360/370  series.  All  complex  variables  are  written 
as  arrays  of  length  two  and  all  complex  arrays  are  written  as  2-dimensional 
arrays  whose  first  dimension  is  two.  This  leads  to  the  code  being  almost 
as  easy  to  read  as  if  the  complex  facilities  were  used.  Complex  division 
is  more  complicated  than  can  be  written  conveniently  as  in-line  code  so  we 
have  taken  this  out  of  line  to  the  subroutines  PA06ED/7ED. 

We  begin  by  describing  the  versions  for  complex  polynomials,  that 
is  PA06AD/8D/CD/DD/ED. 

The  main  call  is  to  a  short  routine  PA06AD  which  calls  PA06BD  to  find 
the  roots  and  PA06CD  to  find  error  bounds.  PA06DD  and  PA06ED  are  short 
subroutines  for  polynomial  evaluation  and  complex  division,  respectively. 

The  time  taken  to  find  error  bounds  is  usually  about  25%  of  that  taken  to 
find  the  roots  themselves  but  can  rise  to  as  much  as  100%  when  there  are 
many  multiple  roots.  Because  of  this  overhead  we  felt  it  was  desirable 
for  the  user  to  be  allowed  to  call  PA06BD  directly  and  instructions  about  how 
to  do  this  are  included  in  the  specification  sheet.  When  calling  PA06CD 
(instruction  9  of  PA06AD)  the  work-space  provided  by  the  user  is  divided 
into  the  four  areas  required  by  PA06CD.  Also  PA06AD  checks  for  zero 
leading  coefficients  and  sets  dummy  error  bounds  to  correspond,  since 
PA06CD  assumes  that  the  leading  coefficient  is  non-zero. 


-16- 


1 


The  algorithm  used  by  PA06BD  to  find  the  roots  has  already  been 
explained  in  section  2  with  the  minor  exceptions  of  the  inclusion  of 


# 


tests  for  zero  leading  coefficients  (roots  at  infinity),  tests  for  zero 
trailing  coefficients  (roots  at  zero)  and  the  scaling  of  all  deflated 
polynomials  so  that  the  largest  coefficient  has  modulus  approximately  equal 
to  the  reciprocal  of  the  modulus  of  the  smallest  non-zero  coefficient. 

We  implicitly  assumed  in  section  2  that  leading  and  trailing  coefficients 
were  non-zero  and  the  inclusion  of  scaling  minimizes  the  likelihood  of 
underflow  or  overflow.  To  avoid  additional  roundoff  we  scale  by  a  power 
of  the  floating-point  base.  For  speed  of  execution  (in  function  calls 
rather  than  basic  arithmetic)  we  use  single-length  working  for  scaling  and 
for  finding  the  initial  iterate.  Also  we  avoid  time-wasting  evaluations 
of  the  moduli  of  complex  numbers  by  working  with  the  squares  of  the 
moduli  or  the  sums  of  the  absolute  values  of  the  real  and  imaginary  parts. 
Subroutines  PA06DD  and  PA06ED,  called  from  several  places  in  PA06BD,  should 
be  regarded  as  part  of  PA06BD.  PA06DD  evaluates  a  complex  polynomial  at  a 
complex  point  and  finds  the  square  of  the  modulus  of  the  result.  PA06ED 
is  a  simple  subroutine  for  complex  division.  The  arguments  of  PA06DD  and 
PA06ED  are  explained  in  comments  at  the  heads  of  each  subroutine.  Other 
coding  details  of  PA06BD  are  explained  in  consents. 

Subroutine  PA06CD,  which  finds  the  error  bounds  using  the  algorithm 
of  section  3,  again  makes  use  of  single-length  arithmetic  whenever  possible. 
In  particular  we  found  that  CABS  executes  significantly  faster  than  its 
double-precision  equivalent  (which  in  any  case  is  not  standard  FORTRAN). 
Therefore  the  error  polynomial  is  held  in  single-length,  the  function 
k(r)  given  by  equation  (3.9)  is  calculated  in  single-length  and  all  the 
Rouche  tests  are  performed  in  single-length.  The  details  of  the  code  are 
explained  in  comments. 
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We  have  tried  to  make  the  code  for  the  case  with  real  coefficients. 


namely  subroutines  PA07AD-PA07ED,  resemble  that  for  the  complex  case  as 
much  as  possible,  using  the  same  labels  and  the  same  comments  wherever  this 
is  appropriate.  The  only  significant  change  in  PA07AD  lies  in  the  way 
the  work-array  W  is  subdivided  when  PA07CD  is  called.  The  only 
significant  changes  in  PA07BD  lie  in  the  code  for  deflation  (instructions 
138-166),  which  is  much  more  complicated  because  we  need  to  test  for  a 
real  root  and  need  to  include  code  for  deflation  with  a  pair  of  complex 
conjugate  roots.  The  code  for  polynomial  evaluation  (in  PA07DD)  is 
longer  because  we  make  use  of  the  fact  that  the  coefficients  are  real  and 
treat  separately  the  case  where  the  point  of  evaluation  is  real.  The 
main  changes  in  PA07CD  are 

(i)  Recognising  complex  roots  when  forming  the  polynomial  n(z-a^.) 
and  multiplying  them  in  as  a  complex  conjugate  pair  to  preserve  real 
coefficients  (instructions  18-28).  It  is  assumed  that  conjugate  pairs 
are  adjacent  in  array  R. 

(ii)  Much  more  complicated  code  (instructions  39-59)  for  finding 
distances  from  the  root  to  all  the  rest.  This  code  avoids  calling 
CABS  to  find  the  distances  between  two  real  roots  and  so  is  much 
faster  where  most  roots  are  real . 

(iii)  Using  the  same  error  bound  for  a  complex  root  and  its 
conjugate  (instructions  114-115). 

6.  Test  results  and  comparisons  with  other  methods 

Our  original  reason  for  providing  a  new  routine  for  the  Harwell 
library  was  that  the  existing  routine  PAOl  sometimes  gave  incorrect  answers. 
The  new  routine  (PA07)  is  able  to  get  answers  all  of  which  have  good  accuracy 
in  about  half  the  time.  We  did  not  pursue  the  error  in  PAOl. 
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We  have  compared  our  algorithm  with  that  of  Jenkins  and  Traub 

(1970)  in  two  ways.  First  we  compared  the  Algol-W  code  given  by  Madsen 

(1973)  with  an  Algol-W  version  of  the  code  of  Jenkins  (1969)  and  found 

Madsen's  to  be  4  to  5  times  faster.  The  required  changes  from  Algol -60 

to  Algol-W  are  very  minor.  Next  we  compared  our  Fortran  code  with  a  rather 

free  translation  into  Fortran  of  Jenkins* code.  We  did  not  try  to  make 

a  literal  translation  but  rather  to  use  his  ideas  to  produce  efficient 

Fortran  code.  The  resulting  program  executed  between  2  and  4  times  slower 

than  ours.  A  further  advantage  of  our  algorithm  is  that  it  is  simpler  and 

2/ 

so  the  code  is  less  bulky  (the  object  code  being  about  3  as  long).  The 
accuracy  of  the  roots  produced  by  the  two  algorithms  was  very  comparable 
except  in  the  last  test  shown  in  Table  1  (Jenkins'  6th  example)  where  we 
obtained  much  smaller  errors  (all  less  than  2xl0’^^)  because  the  roots  of 
modulus  0.9  were  not  all  found  before  those  of  modulus  1  and  therefore  ill- 
conditioned  deflated  polynomials  were  not  generated.  To  compare  the 
error  bounds  with  those  of  Jenkins  we  ran  the  examples  documented  in  his 
report.  In  his  6th  example  our  bounds  were  much  better  simply  because  the 
roots  had  been  found  so  much  more  accurately.  In  his  5th  example  (a 
complex  polynomial  of  degree  21  having  roots  of  multiplicities  1,2  and  3) 
several  of  his  bounds  were  quite  unrealistic  and  our  bounds  were  all  better, 
most  of  them  by  factors  over  100  and  three  (or  seven  if  multiplicities  are 
included)  by  factors  over  1,000.  In  the  remaining  examples  the 
differences  between  the  bounds  were  not  severe. 

Special  purpose  subroutines  exist  in  the  Harwell  library  for  solving 
real  cubic  and  quadratic  equations  (PA03A/AD  and  PA05A/AD,  respectively). 
They  use  direct  methods  involving  only  the  extraction  of  square  and  cube 
roots  but  sometimes  they  lose  accuracy  through  unnecessary  cancellation. 

We  had  hoped  that  a  direct  call  of  PA07BD  might  be  nearly  as  fast  so  that 


we  could  withdraw  the  special  purpose  subroutines,  but  the  speed  difference 
is  by  a  factor  of  about  7.  It  is  hoped  that  better  versions  of  PA03A/AD 
and  PA05A/AD  will  be  written  for  the  library  in  due  course.  Extra  tests 
for  ensuring  reliability  are  likely  to  slow  down  the  programs  a  little,  but 
it  seems  likely  that  good  programs  significantly  faster  than  PA07BD  can  be 
written  for  these  special  cases. 

To  test  our  subroutines  we  read  in  sets  of  roots  and  a  relative  error 
level.  We  used  extended  precision  arithmetic  to  construct  the  corresponding 
polynomial  and  then  made  pseudo-random  perturbations  to  its  coefficients  at 
the  required  relative  level.  The  polynomial  and  its  errors  were  then 
handed  to  our  subroutines.  This  enables  us  to  check  the  actual  errors 

against  the  computed  error  discs.  We  began  by  using  all  the  polynomials  of 

Jenkins  (1969)  and  about  as  many  others  from  local  sources,  but  eventually 
reduced  our  test  set  to  those  shown  in  Table  1  (and  a  few  trivial  ones 
designed  to  explore  corners  in  our  code)  because  the  remainder  showed  no 
useful  additional  information.  We  used  the  same  test  data  for  the  real 

case  by  adding  an  additional  root  consisting  of  the  complex  conjugate  of  any 

complex  root. 

Our  test  results  are  summarised  in  Table  1.  For  each  example  we  show 
the  errors  and  error  bounds  obtained  for  the  first  and  last  root  found 
and  one  intermediate  one.  This  gives  the  reader  an  indication  of  our 
success  in  finding  the  roots  in  order  and  we  have  been  able  to  choose  the 
intermediate  root  displayed  so  that  the  three  roots  together  indicate  the 
full  range  of  conditioning.  We  also  show  the  times  (370/168  secs)  for  a 
call  of  PA06BD/7BD  (roots  only)  and  PA06AD/7AD  (roots  and  bounds).  It  can 
be  seen  that  finding  the  bounds  usually  involves  quite  a  small  overhead. 

The  last  two  cases  are  exceptional  because  of  the  large  number  of  roots  which 
required  discs  containing  many  other  roots.  Such  cases  are  relatively  slow 
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because  we  try  to  find  a  disc  with  only  one  root,  then  try  with  two,  and 


so  on. 

It  was  the  complex  version  of  the  last  example  which  led  us  to 

abandon  finding  disjoint  discs,  because  some  of  the  roots  (e.g.  (0,0,9)) 

are  so  ill-conditioned  that  the  best  disc  obtainable  covers  all  the  roots. 

Therefore  if  we  insist  on  distinct  discs  all  that  can  be  obtained  is  one 

disc  for  all  the  roots.  By  using  overlapping  discs, however,  we  obtained 

12  good  error  bounds  of  which  two  are  displayed  in  Table  1. 

The  real  versions  of  the  third  and  the  last  examples  led  us  to 

abandon  Peters'  and  Wilkinsons'  running  error  analysis  on  the  computed 

polynomial  P(z)  =  n(z-a^).  In  the  last  example  all  the  roots  are  well 

conditioned  but  the  bounds  for  the  errors  in  computing  P{z)  were  so 

pessimistic  that  each  disc  contained  all  the  roots.  Example  3  was  less 

"“8 

dramatic  but  all  the  error  bounds  came  to  about  10  instead  of  about 
10'^^. 
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7.  Appendix.  Specification  sheets  and  listings 

Harwell  Subroutine  Library  PA06AD/BD 

1 .  Purpose 

To  find  all  the  roots  of  a  complex  polynomial 

a-j  +  a2X+. .  .+aj^^^x” 

and  error  bounds  for  these  roots. 

2.  Argument  List 

CALL  PA06AD(A,N.R,E,W,S,NP1,LW)  (double-length  standard) 

CALL  PA06AD(A,N,R,E,W,NPl,W)  (single-length  standard) 

CALL  PA06AD(A,N,R,E,W)  (IBM) 

A  is  a  DOUBLE  PRECISION  (REAL  for  the  single-length  version)  array 
of  dimensions  (2,n+l)  which  must  be  set  by  the  user  so  that  the 
real  and  imaginary  parts  of  a^  are  held  in  A(l,i),  A(2,i).  It 
is  unaltered  by  the  subroutine.  For  the  IBM  version  A  may  be  a 
COMPLEX*! 6  array  of  length  n+1. 

N  is  an  INTEGER  variable  containing  the  degree  n  of  the  polynomial. 

Its  value  must  be  positive. 

R  is  a  DOUBLE  PRECISION  (REAL  for  the  single-length  version)  array 
of  dimensions  (2,n)  used  to  return  the  roots.  These  are  held  with 
real  and  imaginary  parts  in  R(l,i),  R(2,i),  i=l  ,2,...',n.  The  dummy 
value  (1D70,1D70)  is  returned  for  each  infinite  root  (corresponding 
to  a  zero  leading  coefficient).  For  the  IBM  version  R  may  be  a 
C0MPLEX*16  array  of  length  n. 

E.  is  a  REAL  array  of  dimension  at  least  (n+l )  which  must  be  set  by 
the  user  to  error  bounds  on  the  coefficients,  or  to  zero  if  these 
are  accurate  to  machine  precision.  On  exit  the  first  n  locations 
contain  approximate  bounds  on  the  moduli  of  the  errors  in  the  roots. 

W  is  a  DOUBLE  PRECISION  (REAL  for  the  single-length  version)  work 
array  of  dimensions  (2,LW). 
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S  (double-length  standard  version  only)  is  a  REAL  work  array  of 
dimensions  (4,LW),  which  may  be  equivalenced  with  W. 

NPl  (standard  versions  only)  is  an  INTEGER  which  must  be  set  to  n+1 . 

LW  (standard  versions  only)  is  an  INTEGER  which  must  be  set  to 

at  least  5n/4+2  (or  3n/2+2  for  the  single-length  version). 

3.  Alternative  Entry 

The  error  analysis  part  of  a  call  to  PA06AD  takes  typically  about 
20%  of  its  time.  If  speed  is  important  and  error  bounds  are  not  wanted 
then  a  call  of  the  form 

CALL  PA06BD(A,N,R,W,NP1 )  (standard  versions) 

CALL  PA06BD(A,N.R.W)  (IBM  version) 

should  be  made.  The  arguments  are  the  same  as  those  of  the  main  call, 

but  W  need  have  length  only  (2,n+l). 

4.  Method 

The  roots  are  found  by  the  algorithm  of  Madsen  (BIT(1973)  22*  71-75), 
the  principal  features  of  which  are  Newton  iteration  followed  by 
deflation.  The  error  bounds  are  found  by  the  application  of  Rouche's 
theorem  as  recommended  by  Wilkinson  (J. Inst. Maths  Applies. (1971 )  8, 

16-35)  except  that  discs  are  always  taken  with  centres  on  the  approximate 
roots  and  errors  in  multiplying  out  the  polynomial  n(x-R(I))  are  ignored. 
The  disc  for  each  root  is  such  that  it  contains  exactly  the  same  number 
of  approximate  roots  R(I)  as  exact  roots  of  the  true  polynomial.  Note 
that  in  the  case  of  true  multiple  roots  the  corresponding  approximate 
roots  may  be  quite  v/ell  separated  but  each  will  lie  in  the  disc  of  all 
the  others  and  their  mean  will  be  a  good  estimate  of  the  true  multiple 
root. 
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PA07AD/BD 


1 .  Purpose 

To  find  all  the  roots  of  a  real  polynomial 

a^+a2X+...+a^^^x” 
and  error  bounds  for  these  roots. 

2.  Argument  List 

CALL  PA07AD(A,N,R,E,W,S,NP1 ,LW)  (double-length  standard) 

CALL  PA07AD(A,N,R,E,W,NP1 ,LW)  (single-length  standard) 

CALL  PA07AD(A.N,R,E,W)  (IBM) 

A  is  a  DOUBLE  PRECISION  (REAL  for  the  single-length  version)  array  of 
length  at  least  (n+1)  which  must  be  set  by  the  user. to  contain  the 
coefficients  and  is  unaltered  by  the  subroutine. 

N  is  an  INTEGER  variable  containing  the  degree  n  of  the  polynomial. 

Its  value  must  be  positive. 

R  is  a  DOUBLE  PRECISION  (REAL  for  the  single-length  version)  array 
of  dimensions  (2,n)  used  to  return  the  roots.  These  are  held  with 
real  and  imaginary  parts  in  R(1 ,i ),R(2,i ),i=l ,2,. . . ,n.  The  dummy 
value  (1D70,0D0)  is  returned  for  each  infinite  root  (corresponding 
to  a  zero  leading  coefficient).  For  the  IBM  version  R  may  be  a 

t 

C0MPLEX*16  array  of  length  n. 

E  is  a  REAL  array  of  dimension  at  least  (n+1)  which  must  be  set  by 
the  user  to  error  bounds  on  the  coefficients,  or  to  zero  if  these 
are  accurate  to  machine  precision.  On  exit  the  first  n  locations 
contain  approximate  bounds  on  the  moduli  of  the  errors  in  the  roots. 
W  is  a  DOUBLE  PRECISION  (REAL  for  the  single-length  version)  work 
array  of  length  LW. 

S  (double-length  standard  version  only)  is  a  REAL  work  array  of 
dimensions  (2,LW),  which  may  be  equivalenced  with  W. 
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NPl  (standard  versions  only)  is  an  INTEGER  which  must  be  set  to  n+1 . 

LW  (standard  versions  only)  is  an  INTEGER  which  must  be  set  to 

3n/2+2  (or  3n+2  for  the  single-length  version). 

3.  Alternative  Entry 

The  error  analysis  part  of  a  call  to  PA07AD  takes  typically  about  20% 
of  its  time.  If  speed  is  important  and  error  bounds  are  not  wanted  then 
a  call  of  the  form 

CALL  PA07BD(A,N,R,W,NP1 )  (standard  versions) 

CALL  PA07BD(A.N,R.W)  (IBM  version) 

should  be  made.  The  arguments  are  the  same  as  those  of  the  main  call, 
but  W  need  have  length  only  n+1. 

4.  Method 

The  roots  are  found  by  the  algorithm  of  Madsen  (BIT(1973)  13,71-75), 
the  principal  features  of  which  are  Newton  iteration  followed  by 
deflation.  The  error  bounds  are  found  by  the  application  of  Rouche's 
theorem  as  recommended  by  Wilkinson  (J. Inst. Maths  Applies. (1971 )  8, 

16-35)  except  that  discs  are  always  taken  with  centres  on  the  approximate 
roots  and  errors  in  multiplying  out  the  polynomial  n(x-R(I))  are  ignored. 
The  disc  for  each  root  is  such  that  it  contains  exactly  the  same  number 
of  approximate  roots  R(I)  as  exact  roots  of  the  true  polynomial.  Note 
that  in  the  case  of  true  multiple  roots  the  corresponding  approximate 
roots  may  be  quite  well  separated  but  each  will  lie  in  the  disc  of  all 
the  others  and  their  mean  will  be  a  good  estimate  of  the  true  multiple 


SUBROUTINE  PA06A0/A t N,R, t,W ,  NP1«LW) 
SUBROUTINE  PA06AD(A,N«R« E»W ) 

SUBROUrr  INE  PA06A0(AtN(R,EtW.StNPl  tLW) 
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